Terahertz waves dynamic diffusion in 3D printed structures

Applications of metamaterials in the realization of efficient devices in the terahertz band have recently been considered to achieve wave deflection, focusing, amplitude manipulation and dynamical modulation. Terahertz metamaterials offer practical advantages since their structures have typical sizes of hundreds microns and are within the reach of current three-dimensional (3D) printing technologies. Here, we propose terahertz photonic structures composed of dielectric rods layers made of acrylonitrile styrene acrylate realized by low-cost, rapid, and versatile fused deposition modeling 3D-printing. Terahertz time-domain spectroscopy is employed for the experimental study of their spectral and dynamic response. Measured spectra are interpreted by using simulations performed by an analytical exact solution of the Maxwell equations for a general incidence geometry, by a field expansion as a sum over reciprocal lattice vectors. Results show that the structures possess specific spectral forbidden bands of the incident THz radiation depending on their optical and geometrical parameters. We also find evidence of disorder in the 3D printed structure resulting in the closure of the forbidden bands at frequencies above 0.3 THz. The size disorder of the structures is quantified by studying the dynamics diffusion of THz pulses as a function of the numbers of layers of dielectric rods. Comparison with simulations of light diffusion in photonic crystals with increasing disorder allows estimating the size distributions of elements. By using a Mean Squared Displacement model, from the broadening of the pulses’ widths it is also possible to estimate the diffusion coefficient of the terahertz radiation in the photonic structures.


Generalized eigenvalue problem.
In this section we describe our approach to solve Maxwell equations in a patterned layer, specifically a single layer of rods array. For a general incidence condition, the electric field components (α = x, y, z), in the j-th grating layer, expanded in plane waves, read as: E jα (ρ, z, ω) = e ikz ∑ G E jα (⃗ q p + Gx, k, ω)e i(⃗ q p +Gx)ρ (S1) where ρ = (x, y, 0) and ⃗ q p = q xx + q yŷ is the in-plane wavevector. In the following we omit the j index since, in the structures considered, the grating layers are all the same. By magnetic field elimination, Maxwell equations can be written as: and, once Fourier transformed, in mixed coordinates, they turn into: where: q x (G) = q x + G. Moreover, given the dielectric displacement field in the form: The former system can be solved with respect to the z component of the field: The substitution of the expression in Eq.(S6) in the first two equations of the system of Eq.(S5) turns the Maxwell equations in a generalized eigenvalue problem: where the matrix blocks, of dimension (2N + 1) × (2N + 1), are given by : Notice that with the superscript T we indicate vector transposition. In order to transform the system in Eq.(S7) into a canonical problem, we solve the eigenvalue system: Aα n = λ n α n , with: Given U the unitary transformation that diagonalizes A, namely U T AU = λ , by multiplying Eq.(S7) by U T on the left, we obtain: Finally, definingη = λ 1/2η andB = λ −1/2B λ −1/2 we manage to reduce the generalized eigenvalue problem (Eq.S7) to a canonical form: Eventually, the electric-field components (E x , E z ) are computed by transforming the vectorη back to the vector η.

Multilayer solution
In this section we describe how the multi-layer solution can be obtained starting with the simple case of a system composed of two grating layers, as shown in Fig.S1. Moreover we take the in-plane wave vector along the x axis, i.e. we consider the case 2/8 Figure S1. Schematic of the two gratings (j=1,2) configuration with indicated the forward A i and backward B i propagating field amplitudes in the different layers.
φ o = 0. For these conditions the in-plane electric field components in the different regions of the structure can be written as: where the k n 2N + 1 eigenvalues and the η n (G) correspondent 2N + 1 eigenfunctions with 2N+1 G-components are computed for a chosen q x (0) value in the first Brillouin zone (−π/d x , π/d x ).
A phase matrix where ℓ z is the distance between the gratings, matches the fields in the regions of the structure indicated as U and L in Fig.S1: By considering the TE polarization and imposing the Maxwell boundary conditions at z = ℓ z /2 and z = L z + ℓ z /2 we obtain the amplitudes in the transfer matrix form: where: are the interface reflection and transmission amplitudes of the j = 2 grating layer.

3/8
It is, however, well known that 1-3 the inevitable existence of evanescent solutions in dielectric gratings makes the transfermatrix calculation numerically unstable very quickly, so we rewrite the system in Eq.S18 in the scattering matrix form as: η n (G) r > 2 (G, n)A n2 e ik n L z + B n2 e −ik n L z (S18) for the input (A U (G), B 3 (G)) and output (A 3 (G), B U (G)) field amplitudes of the j = 2 grating layer. The forward optical response of the dielectric grating layer, for a given incident in-plane wave vector q p (G ′ ), by the assumption A U (G) → δ G,G ′ and B 3 (G) → 0 makes the system of Eq.S18 a heterogeneous algebraic system that can be solved with respect to the values of the internal electric field amplitudes A n2 , B n2 . Then the system of Eq.S19 gives the matrices of forward reflection B U (G) → r 2 (G, G ′ ) and transmission A 3 (G) → t 2 (G, G ′ ).
Due to the optical symmetry of the grating layer, the backward optical response gives the same reflection and transmission matrices. We can then write the scattering matrix of the system as: Analogous relations can be obtained for the first grating layer: Moreover, Eq.S14 now is: where the field amplitudes are given in the {G, G ′ } reciprocal space, and χ(±ℓ z ) are diagonal matrices. From Eqs.S20, S21 and S25 we obtain: where: The total scattering matrix of the cavity is then: The whole procedure can be iterated when additional grating layers are stacked to obtain the multilayer structure.

ABS optical parameters.
The same experimental set-up used to measure the spectral and dynamical properties of THz pulses transmitted through the photonic crystals was also used to measure the optical properties of the ABS material as a function of frequency. To this purpose transmittance of 3D printed unstructured ABS slabs with thickness of 0.5 and 2 mm were measured. The optical parameters were obtained from the measured spectra by using a numerical procedure able to remove the residual Fabry-Perot oscillations present in the spectral behaviors of samples with flat and parallel surface 4

THz deconvolution
The deconvolution technique was applied to the THz time-domain data both for the ABS slab with a thickness of about 0.41 mm and for the photonic structures made of different layers of grids. In order to recover the impulse response function IRF(t), as a function of the time t, we follow the procedure described in Ref. 6. Results are shown in Fig. S3 and Fig. S4. As a function of the number of layers, the number of peaks in the IRF(t) and their time delay increases (Fig. S4). In addition, an overall increasing of the peaks' width is observed. Since the peaks are overall superimposed on each other, we estimated the width of the first peak appearing from 0.83 ps (1 layer) to 2.3 ps (12 layers), which is quite well separated from other peaks. The results were normalized to the width of the IRF(t) of the THz pulses transmitted in the ABS slab, which reasonably is the width of the incident THz signal. Results are shown in Fig. S5. The behaviour of the normalized pulse width can be well fitted by a linear function 1 + 0.186n ℓ , very similar to that obtained by using the squared module of transmitted THz pulses (1 + 0.193n ℓ ).